use "${clean}newspaper_analysis.dta", clear

// recoding Alexandria, VA to the South (technically part of DC until 1844)
replace south = 1 if city == "Alexandria"
replace north = 0 if city == "Alexandria"
replace state = "Virginia" if city == "Alexandria"

// dropping DC as a border area
drop if state == "District Of Columbia"

// keeping only the north
keep if north == 1

// standardizing variables
egen state_surban = std(state_urban_percent)
egen state_srugged = std(state_ruggedness)
egen state_spostoffice = std(state_postoffice_area)

// keeping the sample consistent
cap drop sample
gen sample = 1
foreach var of varlist $state_npcontrols {
	replace sample = . if `var' == .
}

// creating some variables
local c = 1860

gen period = year > `c'
replace year = year - `c'


// matrices to store the results for graphing
matrix define I = J(6, 3, 0)
matrix define S = J(6, 3, 0)



* regressions
* ------------

xtset pub_id

// lincoln baseline re model
xtreg singular c.year##i.period##ib1.lincoln_county ${state_npcontrols} if sample == 1, re cluster(pub_id)

qui lincom 1.period // intercept change at cutpoint
matrix I[1,1] = `r(estimate)'
matrix I[1,2] = `r(estimate)' + 1.96*`r(se)'
matrix I[1,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix S[1,1] = `r(estimate)'
matrix S[1,2] = `r(estimate)' + 1.96*`r(se)'
matrix S[1,3] = `r(estimate)' - 1.96*`r(se)'
	
local beta_l_re = `r(estimate)'
local beta_l_re : display %4.3f `beta_l_re'
local se_l_re = `r(se)'
local se_l_re : display %4.3f `se_l_re'



// lincoln baseline fe model
xtreg singular c.year##i.period##ib1.lincoln_county ${state_npcontrols} if sample == 1, fe cluster(pub_id)

qui lincom 1.period // intercept change at cutpoint
matrix I[2,1] = `r(estimate)'
matrix I[2,2] = `r(estimate)' + 1.96*`r(se)'
matrix I[2,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix S[2,1] = `r(estimate)'
matrix S[2,2] = `r(estimate)' + 1.96*`r(se)'
matrix S[2,3] = `r(estimate)' - 1.96*`r(se)'

local beta_l_fe = `r(estimate)'
local beta_l_fe : display %4.3f `beta_l_fe'
local se_l_fe = `r(se)'
local se_l_fe : display %4.3f `se_l_fe'	


// mcclellan baseline re model
xtreg singular c.year##i.period##i.lincoln_county ${state_npcontrols} if sample == 1, re cluster(pub_id)

qui lincom 1.period // intercept change at cutpoint
matrix I[4,1] = `r(estimate)'
matrix I[4,2] = `r(estimate)' + 1.96*`r(se)'
matrix I[4,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix S[4,1] = `r(estimate)'
matrix S[4,2] = `r(estimate)' + 1.96*`r(se)'
matrix S[4,3] = `r(estimate)' - 1.96*`r(se)'

local beta_m_re = `r(estimate)'
local beta_m_re : display %4.3f `beta_m_re'
local se_m_re = `r(se)'
local se_m_re : display %4.3f `se_m_re'


// mcclellan baseline fe model
xtreg singular c.year##i.period##i.lincoln_county ${state_npcontrols} if sample == 1, fe cluster(pub_id)

qui lincom 1.period // intercept change at cutpoint
matrix I[5,1] = `r(estimate)'
matrix I[5,2] = `r(estimate)' + 1.96*`r(se)'
matrix I[5,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix S[5,1] = `r(estimate)'
matrix S[5,2] = `r(estimate)' + 1.96*`r(se)'
matrix S[5,3] = `r(estimate)' - 1.96*`r(se)'



* graphing 
* --------
svmat I
svmat S


// temp variable to index estimates
cap drop temp*
gen tempn = _n if S1 != .

gen templincoln = tempn < 3

gen tempfe = mod(tempn, 3)
recode tempfe (1=0) (2=1) (0=.)




twoway ///
	(scatter S1 tempn if templincoln == 1 & tempfe != .,	msize(large) mlwidth(none) mcolor("${orange}") 	msym(D) yaxis(1 2) ) ///
 	(rcap S2 S3 tempn if templincoln == 1 & tempfe == 0,   lwidth(thick) lcolor("${orange}") 	lpattern(dash) yaxis(1 2)) ///
 	(rcap S2 S3 tempn if templincoln == 1 & tempfe == 1,   lwidth(thick) lcolor("${orange}") 	lpattern(solid) yaxis(1 2)) ///
	, ///
	text(.0035 1 "{&beta} = `beta_l_re'") ///
	text(.002 1 "se = `se_l_re'") ///
	text(.0035 2 "{&beta} = `beta_l_fe'") ///
	text(.002 2 "se = `se_l_fe'") ///	
	yline(0, lcolor(gs2) lpattern(dash) lwidth(thick)) ///
	legend(off) ///
	xlab(1 "Baseline" 2 "Fixed Effects") ///
	xtitle("") 	///
	ylab(none, axis(2)) ///
	ylab(0(.005).02, axis(1) grid glcolor(gs7%50)) ///
	ytitle("") ///
	xscale(range(0.5 2.5)) ///
	yscale(range(-0.001 0.016)) ///
	subtitle("Lincoln Counties", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	name(linc_s, replace) nodraw
	
	
twoway ///
	(scatter S1 tempn if templincoln == 0 & tempfe != .,	msize(large) mlwidth(none) mcolor("${green}") 	msym(D) yaxis(1 2) ) ///
 	(rcap S2 S3 tempn if templincoln == 0 & tempfe == 0,   lwidth(thick) lcolor("${green}") 	lpattern(dash) yaxis(1 2)) ///
 	(rcap S2 S3 tempn if templincoln == 0 & tempfe == 1,   lwidth(thick) lcolor("${green}") 	lpattern(solid) yaxis(1 2)) ///
	, ///
	text(.0125 4 "{&beta} = `beta_m_re'") ///
	text(.011 4 "se = `se_m_re'") ///
	yline(0, lcolor(gs2) lpattern(dash) lwidth(thick)) ///
	legend(off) ///
	xlab(4 "Baseline" 5 "Fixed Effects") ///
	xtitle("") 	///
	ylab(none, axis(1)) ///
	ylab(0(.005).02, axis(2) grid glcolor(gs7%50)) ///
	ytitle("") ///
	xscale(range(3.5 5.5)) ///
	yscale(range(-0.001 0.016)) ///
	subtitle("McClellan Counties", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	name(mac_s, replace) nodraw

graph combine linc_s mac_s, ycommon scale(1.5) l1title("Slope Change after 1860", size(small) margin(r=-1)) graphregion(margin(r=-1 l=5)) imargin(tiny)  xsize(100) ysize(50)
graph export "${output}FigD10_np_lincoln_slopes_statecontrols.pdf", as(pdf) replace

	

	
	
	







